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The conditions of multi-phase equilibrium are solved for generic polydisperse systems. The case 
of multiple polydispersity is treated, where several properties (e.g. size, charge, shape) simultane- 
ously vary from one particle to another. By developing a perturbative expansion in the width of 
the distribution of constituent species, it is possible to calculate the effects of polydispersity alone, 
avoiding difficulties associated with the underlying many-body problem. Explicit formulae are de- 
rived in detail, for the partitioning of species at coexistence and for the shift of phase boundaries 
due to polydispersity. 'Connective fractionation' is quantified, whereby one property (e.g. charge) is 
partitioned between phases due to a driving force on another. To demonstrate the ease of use and 
versatility of the formulae, they are applied to models of a chemically-polydisperse polymer blend, 
and of fluid-fluid coexistence in polydisperse colloid-polymer mixtures. In each case, the regime of 
coexistence is shown to be enlarged by polydispersity. 

I. INTRODUCTION 

Materials whose constituent elements are much larger than atoms or simple molecules are now ubiquitous. Such 
substances include polymer blends, colloidal suspensions and emulsions. In the manufacture of such large components, 
it is impossible to produce a pure strain of truly identical particles. Instead, each colloidal latex is minutely different 
in size from its colleagues, and the number of monomers on each long polymer molecule inevitably varies across some 
distribution. These systems are said to be 'polydisperse', though the word is applied only loosely to the polymeric 
case, since not every molecule is unique. Some intriguing experimental |l}-|3| and simulational data have been 
collected from polydisperse systems. For instance, unusual textures have been found in polydisperse polymer blends 
Q , and coexisting phases of hard spheres have been shown to contain unequal proportions of the various particle sizes, 
in both simulation and experiment fl|Jl0||. It is important, then, to understand the statistical thermodynamics 
of polydisperse mixtures. 

It is well known how to calculate the equilibrium state of a thermodynamic system [ pT| . At fixed temperature, 
an expression must be found for the Helmholtz free energy. It is then straightforward to determine any property of 
the system at equilibrium: its value is that which minimises the free energy. To establish the concentrations of the 
system's various constituents, the minimisation must be performed subject to constraints of global conservation. The 
constitution of a uniform system cannot vary, because the amount of each component is conserved, but the free energy 
can sometimes be lowered by partitioning into coexisting phases of differing compositions. 

The non-trivial part of this procedure is to establish a formula for the Helmholtz free energy. In principle, the 
subsequent minimisation to derive phase equilibria is simple. 

For a polydisperse system the story is quite different. Here, the number of components is thermodynamically large. 
In general, this does not greatly complicate the (already difficult) many-body problem of formulating the free energy. 
However, since that free energy is now a function of infinitely many conserved variables, the minimisation procedure, 
though formally understood ]r^ , ^3[ , becomes intractable. 

In the past, many calculations have been performed to study the effects of polydispersity, using pedagogical model 
systems with a convenient form of Hamiltonian [l^Jl5[| , species distribution or approximate free energy 
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3| , p9 -23|. Also, approximate solutions can be found by making ad hoc guesses about the compositions of coexisting 
. £h ! phases, and then minimising with respect to just a few variables such as the overall density of each phase |TJ,Q. Al- 
ternatively, the minimisation can be performed numerically, by replacing the continuous distribution of concentrations 
S-h ' with an arbitrary, finite set of variables |2f| . 

The above methods invariably involve some way of approximating the many-body problem inherent in the thermo- 
dynamics of the particular system. The present study differs in that wc shall calculate only the effects of polydispersity. 
The analytical procedure was outlined recently |pfj-|28||, and is presented in more detail in the present article, and 
applied in new ways to solve a range of problems in the physics of polydispersity. The method is a perturbative one, 
whereby polydisperse phase equilibria are derived from a reference state that is monodisperse and also in multi-phase 
equilibrium. We make no assumptions as to the nature of the system, and require no special properties or approxima- 
tions to ensure tractability of the statistical mechanics, since we shall not address the underlying many-body problem. 
Instead, we assume that the many-body properties of the monodisperse reference system are known (whether from 
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experiment or theory), and derive some exact thermodynamic relations, applicable to most systems, for the changes 
in its properties due to the introduction of polydispersity. 

The aim then is to treat a polydisperse system as a perturbation to a monodisperse reference state, using the 
variables that characterise the polydisperse particles (for instance, size and charge) as expansion parameters. There 
are two major obstacles to such a theory. Firstly, the reference state is singular. To calculate thermodynamic 
equilibria, a knowledge of the chemical potential of each species is required. In the reference system, the population 
is zero for all species but one. Unfortunately, the chemical potential of a non-populated species is (as a rule) negative 
infinity. The second obstacle is that a narrow distribution of species is not necessarily a smooth one, and therefore 
may not be conducive to small-variable expansion. To proceed, then, we must isolate the badly behaved functions for 
exact treatment, and expand only smoothly varying quantities. The distribution of species is not expanded (contrast 
p6|), so that any such distribution may be treated, including a set of delta peaks representing a finite number of 
components. The canonical ensemble is used throughout, as it is of greatest practical relevance. Thus, the overall 
mix of species (the 'parent distribution') is specified a priori, and the pressure and chemical potentials are derived 
quantities, used to calculate the resulting distribution adopted by each coexisting phase after partitioning. 

Perturbative methods suggested previously [ p9|]3C}l have been dogged by a proliferation of variables, leading to 
unwieldy expressions. That problem is avoided here by relating all properties of the polydisperse system to just 
two or three functions that parameterise a generic free energy. An alternative approach for perturbing about a 
monodisperse reference is to assume that all species but one are dilute, and to use their concentrations as small 
expansion parameters fl3| , in the manner of a virial series. In the present study, by expanding in the width, rather 
than height, of the distribution of polydisperse components, we are using a reference state that correctly treats the 
many-particle interactions of concentrated, phase-separated systems. For narrow distributions therefore, this series is 
expected to converge more quickly than a virial-like expansion. An exception is at critical points, where both methods 
fail. 

The formalism and notation are set out in the next section, where a generic free energy is perturbatively expanded 
to low order, to derive the pressure and chemical potentials of a single polydisperse phase. These expressions are used 
in section [IlJ to analyse polydisperse phase equili bria. The normalised distributions, giving the relative amount of 
each species in each phase, are calculated in section III A . Their differences (due to fractionation) are found to obey a 
very simple law. For systems that are simultaneously polydisperse in two properties (e.g. size and charge), this law is 
shown to describe 'conve ctive fractionation', whereby one property is partitioned between phases due to a driving force 
on the other. In section III B , the effect of polydispersity on the total number densities at coexistence is calculated. 
The resulting shift of a phase boundary along the density-axis of a phase diagram is shown to be proportional to 
the variance (standard deviation squared) of the parent distribution, in the limit of a narrow distribution. Explicit 
expressions are given for the cloud- and shadow-point densities, as well as more general binodals. All the derived 
formulae quantifying polydisperse phase equilibria are expressed in terms of a couple of basic parameters of the free 
energy. In case the free energy is not know n for the polydisperse system in question, a meth od for evaluating the 
relevant parameters is given in section [VA for soft interaction potentials, and in section IV B for hard spheres and 
related systems. The method for soft potentials is extended to determine the lowest-order effect of polydispersity on 
correlation functions in appendix |A|. To demonstrate the ease of use of the formulae for phase equilibria, they are 
applied in section ^ to a Flory-Huggins model of chemically polydisperse polymers, and in section VI to fluid-fluid 
coexistence in colloid-polymer mixtures, where slight polydispersity is shown to widen the coexistence region. 



II. PROPERTIES OF A SINGLE PHASE 



Consider a polydisperse phase of overall number density p, whose constituent particles are distinguished by a number 
v of properties (e.g. mass, charge, diameter, oblateness. . . ). Each particle is characterised by an ^-component vector 
e, drawn from a narrow distribution p(e) normalised over the whole of e-space. Let each component of e represent 
the (dimensionless) fractional deviation of each property from some reference value. For instance, a hard sphere of 
radius r in a polydisperse sample would be characterised by a scalar e = (r — ro)/?"o, in terms of the reference size r$. 
Since the distribution p(e) is narrow, the components of e are always small numbers. 

To isolate the badly behaved part of the chemical potential, we shall write the phase's free energy density (which 
is a functional of the distribution of species densities pp(e)) in two parts 

%( £ )]Er d + r x . (i) 

This is an exact statement, as it serves only to define T evi as the excess free energy density, over and above that of a 
polydisperse ideal gas of densities pp(t). The ideal part is simply the sum, over the continuum of species, of the free 
energy density of an ideal gas of each species, 
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T id [pp(e)} = J d<e pp(e) {Hpp(e)) 1} . (2) 

The chemical potentials for the continuum of species are given by a functional derivative |l2| of T with respect to the 
density of each species, 

M(C) " We)] (3) 

which, from Eq. ([!]), splits into two parts, 

p(e)=p id (e) + p cx (e). (4) 

By functional differentiation of Eq. (Q), the ideal part of the chemical potential is 

M id (e) = Hpp(e)) (5) 

which is singular in the monodisperse limit p(e) —* (5(e), because it derives from the entropy of mixing. The other 
part, pf x (e), remains well behaved, since it describes the physics of interactions which vary little from one species to 
the next. 

The obstacles to a small-variable expansion, discussed above, have now been isolated in the badly-behaved but 
well-defined function p ld (e) in Eq. (|5|), and the narrow, but possibly rapidly varying distribution p(e), also appearing 
in Eq. (||). These quantities will be treated exactly, while the effects of interactions, embodied in the excess free 
energy density T cx , are expanded in the small variable e. 



A. Expansion of the excess free energy 

The excess free energy density in a phase at equilibrium is a function(al) J- cx ([p(e)], p) of the intensive variables 
p(e) and p. The dependence of !F ex on temperature and other fields is suppressed in the notation, and T cx is assumed 
to be absolutely minimised with respect to any non-conserved order parameters. 

Any distribution p(e) is uniquely defined by its moments {(e) , (ee) , (eee) , . . .} which, for vectorial e, are averages 
over p(e) of outer products of e. For instance, 



(ee) = J eep(e)de. 



Since p(e) is narrow, higher moments are increasingly small. So a controlled small-variable expansion of T cx is 

T cx (p, (e) , (ee) ,...) = T cx (p, 0, 0, . . .) + (e) • A(p) + (ee) : B{p) + (e) (e) : C(p) + 0(e 3 ) (6) 

where F ox (p, 0,0,.. .) is the free energy of a phase of monodisperse particles all characterised by e = 0. We shall write 
this as ^f(p), where the subscript m denotes a quantity in the monodisperse reference system. The coefficients A, 
B and C are vector and tensor functions of the overall density p. The notation 0(e 3 ) indicates terms of third order 
in e and higher. Equation (^]) is a generic expansion, not specific to any particular system. The use of a non-singular 
expansion, involving only integer powers of small quantities, has the status of a conjecture. An alternative motivation 



for the form of Eq. (H) is given in section IV 



B. Expansion of the excess chemical potential 

Following Eq. (§, Eq. (|) is differentiated, using the identity 

8T CX _ dT cx dT cx 5 (e) dT cx 8 (ee) 



S[pp(e)} dp d(e)8[pp(e)} (ee) S[pp(e)} "" 

where 5 (e m ) /S[pp(e)] = (e rn — (e m ))/p. In this notation, dT ex /d (e) is a vector whose components are the derivatives 
of T eyL with respect to each component of (e). Similarly, dT cyi jd (ee) is a tensor. Thus we obtain a general expression, 
truncated at second order in e, for the excess chemical potential p ex (e) of the species with property e, in a phase 
characterised by {p,p(e)}, 

^ ox (e) = /C(p) + (e) ■ (A'(p) - A/p) + (ee) : (B'(p) - B/p) + (e) (e) : (C(p) - 2C/p) 

+ e- (A + 2(e) -C)/p+ ee: B /p + 0(e 3 ) (7) 

where p^(p) is the excess chemical potential, dT^ /dp of a monodisperse reference phase at density p. Equation (0) 



will be central to the calculations of phase equilibria in section III 
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C. Expansion of the pressure 



The pressure of a polydisperse phase with density p and distribution p(e) can be found from 

P = -T+ J n{e)pp{e) d u e 

which is the continuum analogue of the standard thermodynamic relation for the pressure of a mixture. Using Eqs. (|l|) 
and (f|) it follows that 

P = P id -T cx + pJ p cx (e) p(e) d v e (8) 

where P ld is the ideal pressure of the polydisperse mixture. It is easily confirmed from Eqs. (||) and (||) that this 
is the usual ideal gas pressure P ld = p. Substituting the series expressions for the excess free energy and chemical 
potential (Eqs. ^ and ^7§) into Eq. (||), and using the pressure of the monodisperse reference system at density p, 
Pm(p) = P — 3~m(p) + pfi(p), yields the expansion for the pressure of the polydisperse phase, 

P = P m (p) + (e) ■ {pA'(p) -A) + (ee) : ( P B\p) - B) + (e) (e) : (pC(p) - C) + 0(e 3 ). (9) 



III. PHASE EQUILIBRIA 
A. Fractionation 

At coexistence between two or more phases, the distribution p{e) of properties is typically different for each phase. 
We shall assume that the overall distribution of all particles in all phases is known (since the particles might have 
been synthesised en masse, before phase separation was induced). This will be called the parent distribution p p (e). 
Without loss of generality, the origin of e will henceforth be set at the mean of this parent distribution, so that 
(e) = 0. The aim of this section is to calculate the distribution p a (e) in any given phase a, of M coexisting phases. 

The number of particles in phase a is some fraction n a of the total number in the system. So conservation of 
material is expressed as np = 1. In fact, a stronger criterion than this holds. The number of particles belonging 
to each species is conserved. That is, 

M 

^npppie) =p p (e) (10) 



for all e. 

This is a convenient point at which to introduce a special notation which will be useful later. If Q is any property 
of a phase then, for that phase, let A[Q] be its deviation from the mean of Q over all phases, weighted by the number 
of particles in each. Hence, in phase a, 

M M 



The notation [. . .]jg indicates the difference between the quantity evaluated in phase a and in phase f3. 

Returning to the problem of coexistence, note that, in any two of the M coexisting phases, a and (3 say, the chemical 
potential is equal, i.e. p a (e) — pp(e) for all e. Combining Eqs. (Q) and (||), this implies 

p a p a (e) exp^ x (e) = pppp{e) exp^f(e) (11) 

which can be used to eliminate each pp(e) in Eq. ([l0|), yielding 

P P ( e ) 



Pa{e) = 



E/ji/^exp p™(e)-pf(e) 



(12) 
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We shall now use Eq. (|) for the excess chemical potential, truncated at first order, 

p cx (e) = vg(p) + <e) ■ (A' - A/p) + e ■ A/ p + 0(e 2 ) 

and expand the exponential in Eq. as 

exp[^»(e)]g = (1 + e • [A/p]£)exp + (e) ■ {A' - A/p)] a p + 0(e 2 ). 

So Eq. ( O ) can be written 

/ \ ft,(e) A gjA^g exp [pg + jj ■ (A' A/p)]; \ 
p a {e) = I 1 - e + 0(e )\ (13) 

where <S Q = X)/3 n /3^ ex P[Mm + ( e ) ' — A/p)] p. If Eq. ( |l3| ) is integrated with respect to e, then the term linear 
in e disappears, due to the condition (e) = 0. We are left with 1 = 1/S a + 0(e 2 ), or S a = 1 + C(e 2 ), which 
can be substituted back into Eq. (|l3|). Finally, integrating Eq. (|Tl| ) with /i cx (e) expanded to zeroth order yields 
{Pa/ Pp) ex P [^m]/3 = 1 + C(e), which we also substitute into Eq. (|13[). The result is 

p a (e) = p p (e) {l - e ■ A Q [A(p)/p] + 0(e 2 )} . (14) 

This is the result we wanted: an expression for the distribution of species in each phase. Not surprisingly, it is 
almost the same as the parent distribution (which remains unapproximated) , since all species have similar properties 
(e is small for each). The lowest-order correction to p P (e) is expressed in terms of A, a (vectorial) parameter of the 
excess free energy (Eq. (^)), and pp, the density of each phase. In principle, the quantities on the right hand side 
of Eq. (|l4|) are known. The densities of the phases p a will differ a little (by an amount 5 a , say) from the binodals 
of the monodisperse reference system p™, so that p a = p™ + 5 a , and the numbers of particles n a will have similar 
small changes. However, in Eq. ( p^ ) one may substitute either the monodisperse or the polydisperse values of p and 
n, whichever are available from theory or experiment, since the small differences affect only higher-order terms. 

A more elegant expression is obtained if we find the difference between the normalised distributions pie) of two of 
the M coexisting phases. From Eq. (fl4|), 

\p(e)]^-p p (e)e-[A/p] a (15) 

which, together with Eq. (|Tc|) , contains the same information as Eq. (fl4|). The arrow(— ►) indicates the strict asymptotic 
limit as (e 2 ) p — * 0. 

The parameter A/p in Eq. (|l5|), is the lowest-order coefficient of e in Eq. (|7|), so we can write 

[p(e)]2--p p (e)e- \v e p™(e) V. (16) 

A gradient in e-space is denoted V^, and it is evaluated for phases a, (3 in Eq. (|l^) at the mean species e = 0. 
Equation ( |T6| ) says that, at coexistence, the distributions separate along the steepest gradient of /Lt„ x — pJ^ in e-space. 

It is informative to take moments of Eq. (|l^), since these will tell us how the various properties of the particles are 
correlated within the phases. One might wish to examine averages of the individual components £j of e, such as (ei) 
the mean value of some property such as the size of particles, or (£162) the cross-correlation between e.g. sizes and 
charges of the particles. The most general moment, from which this information can be extracted, is the mean of an 
outer product of m vectors e. That is the mth rank tensor 

( e m ) = J ££_^_£ p(e) d v e (17) 

771 

whose components are all possible correlators of order m. From Eq. (|l5|), the difference of such a moment between 
two of the coexisting phases is 

KOI; = - (£ m+1 ) p • [A/p]; + o(£ m + 2 ) (is) 

which involves the next moment of the parent distribution, and one summation (the scalar product) over the v 
polydisperse properties. The scalar version of Eq. ( n8[) w as derived previously for substances that are polydisperse in 
a single property, at two-phase pq] and multi-phase |27| coexistence. 
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It has been noted (3^] that for some distributions, odd moments ^e 2fe+1 ) p are not of order (2fc + l) in the width of the 
distribution but, due to cancellation in the integral (Eq. (p7|)), are of higher order, so that the RHS of Eq. ( |l8| ) is zero 
to 0(e m+1 ) and the unevaluated 0(e m+2 ) terms are the dominant ordeqj. Specifically, this is the case for distributions 
which become symmetric in the narrow limit (e.g. Gaussian, Schultz). For strongly asymmetric distributions though, 
excepting unlikely cancellations, the mth moment is of mth order in the width. Equations ( [To] ) and (|l^) remain valid 
for all narrow distributions. 

As the notation of Eq. ( |l8| ) is rather abstract, let us consider a specific example: a system of particles with a 
range of radii r and charges q. At equilibrium the particles are partitioned into two coexisting phases denoted a 
and P, and we are interested in the difference between the average radii in these phases. The overall average radius 
throughout the system is called (r) p and the average charge {q) p . Applying Eq. ( |l8| ) with m = 1, t\ = r / (r) p — 1 
and 62 = qj (q) p — 1, and taking the first component gives the answer 

(r) a (r) = a [(r 2 ) p (r)j] + b [(rq) p (r) p (q) p 

with coefficients a = — [Ai/p]^/ (r) and b = — \A-ij p^j (q) . The first term says that, in the absence of charges, the 
amount of size fractionation is proportional to the overall size variance (the square of the standard deviation), with a 
coefficient a indicating the difference between each phase's affinity for large particles. With charge polydispersity, the 
second term indicates that size fractionation can occur even if a = so that neither phase favours larger particles. If 
size and charge are cross-correlated in the parent, i.e. bigger particles tend to have bigger (or smaller) charges, then a 
'convective' fractionation of sizes will be driven by b, the 'driving force' separating charges. These two contributions 
are additive for a system with a narrow distribution. 



B. Shift of the phase boundaries 

So far we have established the lowest order change (with width of the distribution) of the distribution of properties 
in each phase, relative to the original, 'parent' distribution (Eq. (Q)). However, a knowledge of the normalised 
distributions is not a complete characterisation of the coexisting phases. We would also like to know the overall 
densities of the phases, p a . For a monodisperse system, the binodal densities are independent of the number of 
particles in each phase. This is not so when the system is polydisperse, since the size of one phase's share of the 
available particles influences the shape of its distribution of species and, hence, its thermodynamic properties. One 
coexistence of particular interest is the 'cloud point'. A phase is at its cloud point when its density lies on the 'cloud 
curve' in the phase diagram. It then coexists with an infinitesimal amount of another phase, whose density is on 
the corresponding 'shadow' curve. The cloud point is important from a practical point of view, since it marks the 
threshold at which a separation of phases is first observable, and the resulting morphology and fractionation can have 
useful applications. The cloud point also inspires theoretical interest since it is relatively easy to treat in a two-phase 
system, where the distribution of species in the majority phase is simply equal to the parent. 

In the present analysis, as in the previous section, we shall calculate the conditions of general phase equilibria, not 
just the cloud and shadow curves. Having found the general solution, we shall return to the special case of the cloud 
point. 

To establish the phase equilibria, we shall demand that the pressure and chemical potentials are equal in all phases, 
using the expressions for pressure and chemical potential derived in section |J. Before proceeding further, let us take 
moments of Eq. ([l4|), as this will lead to some substantial simplifications. We see that the second moment, (ee), in 
any given phase is, to lowest order, equal to that in the parent. However, the mean, (e), in that phase is also of second 
order in the parent's width, since the mean of the parent vanishes by definition. In orders of the width of the parent, 

<jp = (e ■ e) p , we have 

(e) = - (ee) p • A [A/p] + 0(a 3 p ) (19a) 
and (ee) = (ee) p + 0(o*). (19b) 



1 When m — 2, Eq. ( |l8| ) relates the second moment of the daughters to the skew of the parent |26|. For highly symmetric 
distributions the skew is small (or vanishing), so the LHS of Eq. ( |l8| ) is small, though not vanishingly so, as it is limited by 
the unevaluated C(e m+2 ) terms. On the other hand, for a strongly asymmetric parent, by Eq. (|l8|), the daughters have very 
different second moments. Thus the skew of the parent is highly influential, contrary to the impression given in Ref. ]3l|. 
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So the terms involving (e) (e) in Eqs. (R) and can be dropped, since they are of order a^. Let us clarify this point. 
The excess chemical potentials in Eq. (Q) and the pressure in Eq. (|^) were calculated consistently to second order in 
the deviations e of the properties of a phase's constituent particles from some reference. In general, then, to evaluate 
these thermodynamic potentials for an arbitrary phase, with respect to an arbitrary reference, terms involving (e) (e) 



must be preserved. However, the choice (e) p = 0, to fix the arbitrary reference, allows us to drop the terms in Eqs. (0) 
and ^ with coefficient C , for coexisting phases. This is because we have found that when several phases coexist, their 
mean properties differ very little from each other, so that (e) Q — (e) p is second-order small. 

To calculate the densities of the coexisting phases, we first demand that the pressure P(p, (e) , (ee)) is equal in any 
two phases a, (3 of the M coexisting phases. We shall write this as 



P(p, (e) , (ee)) 



= 



(20) 



and use the expression for P(p, (e) , (ee)) in Eq. (|). Equation © is a condition on p Q , the binodal densities which 
we want to find. As mentioned above, in the limit of a narrow parent, these densities are close to the monodisperse 



binodals p™, differing by an amount S a , so that p a = p™ + S a . We shall determine the small shift S a . In Eq. we 
can Taylor-expand the monodisperse pressure P m (p a ) — P m (Pa) + i^(p") It will transpire that S a is of order 
<7p, so that it is sufficient to truncate this Taylor expansion at first order in S a . The condition of pressure balance 
(Eq. (p0|)) becomes 

[Pm{pm) + PLiPm) 6 + (e) • {pM - A) + (ee) : {pB 1 - B)]% = 0. 

We see that, due to the low order of expansion, the change in density 6 affects only P m , the 'monodisperse 
contribution' to the pressure. Effectively, each of the polydisperse phases can be replaced by a monodisperse phase 
that is subject to an 'external field' ((e) ■ (pA 1 — A) + (ee) : (pB 1 — B)). In this interpretation, each (fictitious) 
monodisperse phase responds to the field with a density change 5, governed by the response function P^. The same 
interpretation can be applied to the chemical potential balance below. 

We can eliminate the first term because the monodisperse reference phases at densities p™ are in pressure balance, 
that is [P m (p m )]a = 0. The density derivative of the pressure in the monodisperse system, P^(p m ), can be expressed 
in terms of the phase's isothermal compressibility, 



-1 



dV 



N,T 



dP 

dp 



dp 



Hence pressure balance dictates 



— + (e)-(A'p-A) + (ee) : (B' p - B) 
up 



= 



(21) 



(22) 



for any pair (a, (3) of the M phases. 

Another constraint on the shifts in density 5 a .p can be derived from the equality of all chemical potentials between 
any pair of phases, which leads to Eq. ( fll| ) for the excess parts. We substitute from Eq. (R) for p J cx (e), and further 
substitute 

/C0») = MmO) - Inp 

= - lnp + p m (Pm) + Mm(Pm) 5 

to yield 

\p(e) exp {p m (p m ) + p[J + (e) • (A' - A/p) + (ee) : {B 1 - B/p) + e • A/p + ee : B/p}]^ = 
Expanding the exponential of small quantities and integrating gives 

[e^)(l + p[J + (e) • A' + (ee) : {B 1 + A A/2p 2 })]^ = 0. 
Now, applying Eq. ( pl| ) and the condition [p m (Pm)]a = yields 

5 



np^ 



(e) -A' + {ee) : (B' + AA/2p 2 







(23) 
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which, with Eq. (p2|), forms a pair of simultaneous equations in two unknowns, S a and 6p. 

As an aside, it is interesting that, as well as fixing the M unknowns S a , the 2(M — 1) constraints in Eqs. ( p2|) and 
( p3| ) impose extra conditions on the reference system for M > 2. The extra constraints are analogous to the Gibbs 
phase rule: a single component, whether monodisperse or near-monodisperse, can exist as a pair of coexisting phases 
under a range of conditions, whereas triple or higher coexistence requires a special choice of external fields [B2|. It is 
intriguing that the extra conditions imposed on the monodisperse reference system constrain the values of A{p) and 
B(p) (defined in Eq. (^)), which, though defined in the monodisperse system, are not manifestly relevant to its phase 
equilibria. 

It is straightforward to solve the linear equations (E2 and E3I) to yield the value of S a , expressed as 



IvYc 



(24) 



in terms of x = 5/np+ (e) • (A' p - A) + (ee) : (B' p - B) and of y = (e) ■ A/p+ (ee) : (B / p + AA/2p 2 ), where ex- 
pressions for (e), (ee) can be substituted from Eqs. fl9|). However, the expression is inelegant, since there are M 
phases, of which (3 is in no way special. The shift S a in the a-binodal should have a symmetric dependence on all 
other phases. Since (3 is arbitrary, we may symmetrize Eq. ( p4| ) by averaging it over all /?, with respect to any weight 
factor z a p\ 



•"a v-^ 

For neatness, we choose z a p = np (p^ 1 — Pp 1 )- 

The final answer for the shift 5 of the density of oneR of the M coexisting phases due to polydispersity is 



A 



A[A/p]A[A/p] -2A[B/p] 



S = pn <ee) p : ( (A' p - A) A[A/p] B' p + B + -J= ) ( 25 ) 

plus terms of order a^. As this is a lowest-order expression, the densities p appearing on the R.H.S. may be taken as 
either the monodisperse or the polydisperse binodal values, whichever is most convenient. 

As yet, we do not know the values of the parameters A, B, A' and B' in Eq. (p5[). They will be discussed in Section 



[V. Nevertheless, the structure of Eq. ( |25| ) is informative. As stated above, S is indeed of order a 2 . If, as is often the 
case P,[7|, p3| , p4||3^j3^ ] , a phase diagram is drawn for some system, showing the densities of coexisting phases against 
the width cr p of some parent distribution, then Eq. ( |25| ) gives the leading-order Taylor expansion of the shape of the 
phase boundaries^. By definition, they meet the (er p = 0)-axis at the monodisperse values p m , and we have shown 
that they do so perpendicularly, since there is no term linear in a p . According to Eq. ( |2^ ) , the curvature of each phase 
boundary is proportional to the isothermal compressibility k of the phase. Hence the small-variable expansion breaks 
down at a critical point, where the compressibility diverges. This is not unexpected, since the position of the critical 
point can itself be shifted by the presence of polydispersity, for instance to a different temperature. Criticality in the 
monodisperse system is a coexistence between two (or more) phases of equal densities. Under the same conditions, 
slightly-polydisperse phases, being non-critical, will either coexist at finitely-different densities, or not demix at all. In 
either case, the thermodynamic state is far from the monodisperse reference, so the divergence of Eq. ( p5[ ) at criticality 
is correct. 

Equation ( p5| ) exhibits a property typical of polydisperse systems. For a given parent distribution, the positions of 
the phase boundaries depend on the number of particles in each phase (since this affects the averages in the definition 
of A[. . .]). This contrasts with monodisperse systems, for which tie-lines can be drawn, connecting a set of points 
in the phase diagram's forbidden region to the same coexisting end-points. It is interesting that this signature of 



2 The subscript a, denoting the phase in question, is now dropped since it applies to every quantity in Eq. ( pq ) (except (ee) p 
and operands of A). 

3 This expansion of the equations of the phase boundaries is in some ways equivalent (though in a quite different formalism) 
to the Clausius-Clapeyron-like equation derived by Bolhuis and Kofke || for the P — v phase diagram, where v parameterises 
the variance of a Gaussian distribution of activities of polydisperse hard spheres. The present analysis, though restricted to 
narrow distributions, is more general. 
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polydispersity p3] is evident at the lowest order in a p . Hence, as soon as any shift in the phase boundaries due to 
polydispersity is observable in a system, that shift should vary as a tie-line is traversed. 

As mentioned above, a coexisting state of particular interest is that between an infinitesimal amount of one phase 
(the 'shadow' phase) and a majority phase defined to be at its 'cloud' point. Such a state defines the extreme 
boundary of the coexistence region. The movement of the cloud point^, 6 C is of great importance, since this tells an 
experimenter where, in the phase diagram, to first expect phase separation to occur. Let the portion of particles in 
the shadow phase be the infinitcsimally small number n s = s. Then the fraction belonging to the majority, cloud 
phase is n c = 1 — s. For any property Q, A[Q] vanishes in the majority phase as s — > 0, since the cloud phase defines 
the number-weighted mean. In the cloud phase, 

A C [Q] = Q c - (n c Q G + n s Q s ) = s [Q] c s 

where [Q]° denotes the difference between the values of Q in the cloud and shadow phases, Q c — Q s . In the shadow 
phase, 

A 3 [Q] = (l-a)[Q]' c . 
Substituting these expressions into Eq. (p5|) gives the shift of the cloud point as 



6 C = -p c K c (ee) p : [ B c p c - B c + 2[l/p] c ' ' 



The coexisting shadow phase has a density p" 1 + 5 S where 

Ss = -p.«. (ee) p : Ps -B s+ ^AM^^M + (A ; pj _ Ag) [A/p]; ) 



(27) 



It is not always appropriate to express the position of a phase boundary in terms of the number density of particles, 
p = p m + 5. For instance, the phase diagram is often represented in terms of the volume fraction of particles, </>. 
This is the product of the number density and the volume of a particle. Since fractionation occurs, so that the mean 
volume of a particle differs from phase to phase, (f> and p are not equivalent, and the coexistence region has a different 
shape in the two representations of the phase diagram. Quantities such as <f> are easily derived from the above results 
as follows. For spherical particles, for instance, using the scalar e to denote fractional deviations in particle radii: 
r,; = (1 + e,;)r m , we have 

^=^<r 3 p) = ^((l + e) 3 )(p m + 5) 
=> + — +3<6)+3<6 2 ) +0(al) (28) 

(Pm Pm 



with the quantities on the RHS given in Eqs. (|19|) and (25) 



In describing the phase diagram, one final quantity of interest is the width of the coexistence region, sometimes 
called the miscibility gap. In terms of p, this is the range of densities of the system as a whole, for which it separates 
into more than one phase. Say the miscibility gap is bounded by two phases, called a and (3 in the reference system. 
Quenching to one edge of this gap, the cloud point of phase a, will cause an infinitesimal amount of phase (3 to form. 
At the other extreme, a vanishing amount of phase a coexists with majority phase 0. So the miscibility gap is given 
by application of Eq. ( |26| ) alone, to establish the two cloud points of phases a and /3. The miscibility gap [p c ]p is then 
given, in terms of the gap in the monodisperse reference system, by 

[Poll = [P?] a + : { M/J [B/P\ a p + [l/p]% \P<B' P B))1 + [A/p]% [A/ p]% ( Pa n a + Pf} K )/2). (29) 

Note that, with the gap defined to be positive, [— is positive. In principal, the lowest-order change in the 

miscibility gap due to polydispersity, given in Eq. (gg), may be positive or negative, depending on the system. We 
shall return to this topic in section [y]. 



4 Greek subscripts were used to denote generic phases, whereas c, s indicate cloud and shadow phases, for which the relative 
numbers of constituent particles are specified. 
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IV. EVALUATING THE PARAMETERS 



We have determined the state of thermodynamic equilibrium for polydisperse phases, in terms of certain properties of 
the monodisperse reference phases. Some of these properties, the density p and compressibility k, are well established 
for most common systems. Others, the parameters A(p), A'(p), B(p) and B'(p), require further discussion. These 
parameters appear in the excess free energy, Eq. (|^), and we require their values in the monodisperse system (where 
(e) = (ee) = 0) at some density p. One might already have an expression to hand for the free energy of a particular 
system, in the form of Eq. (^), in which case it is trivial to identify the appropriate parameters and substitute them 
into the various expressions above. The results for two commonly-studied systems are given in sections and [v| 
Alternatively, the required numbers may be established by experiment, or may require derivation from first principles. 

We shall now find some generic expressions for A(p), B(p) and, f or completeness, also C(p), which will inform 
our physical interpretation of these quantities 



[V A| , we shall apply those expressions, to construct 

That theory 



Then, in section 

a method for first-principles calculation of the parameters from thermodynamic perturbation theory, 
is applicable whenever the Hamiltonian of the system can be differentiated with respect to the properties e of the 
particles. Though many systems meet this criterion, one system of particular theoretical and pedagogical interest, a 
set of hard spheres, does not. Its Hamiltonian is discontinuous, being zero for all physical configurations and rising 
discontinuously to infinity if any spheres overlap. However, hard spheres belong to a special class of systems (scalable 
systems) for which the parameter A, which co ntrol s the emergence of fractionation (see Eqs. (|l4|, |l5, 18)) is calculable 



by an alternative method, detailed in section IV B 



To find general expressions for the quantities A, B and C, let us Taylor expand a generic excess free energy with 
respect to the properties 6j of each constituent particle separately, where i labels the N particles. We expand to 
second order in these N vectorial variables thus 



F c 



N 



dF c 



N 

e ^ : 



d 2 F c 



06; 9e, 



0(e 3 ) 



where the subscript m indicates evaluation in the system where e, = for all i. In that system the particles are all 
the same, so the derivative dF cyL /dei takes the same vectorial value for any particle i. Hence we may evaluate it for 
particle number 1, without loss of generality. Similarly, the second derivate d 2 F ex / deidej takes just two tensorial 
values, depending on whether i and j label the same particle (number 1, say) or different particles (1 and 2, say). 
Hence, dividing the above equation by the system volume V gives the excess free energy density 



dF° 



1 . . d 2 F CK 

nP\ ee ) : q — 5 — 

2 oe\de\ 



1 d 2 F CK 
+ -p(e) (e):N—— 

2 oeide2 



0(e 3 



(30) 



Note that, although the final term contains a factor N, it remains intensive, because d 2 F cx /deide2 ~ 1/V since 
particles 1 and 2 interact less as the system grows. 

Comparing Eq. (^) with Eq. (@), we can identify expressions for the coefficients A, B and C, which we also express 
in terms of derivatives of Eqs. (g) and ([?]), 



C = 



dT cx 

g-pex 



d{ee) 
1 d 2 T ex 
2d(e)d(e) 



A = 
B = 



= P- 



dF c> 
dei 



1 d 2 F c 



= P- 



dp cx (e) 



= nP 



1 dV x (e) 



2 deidei 



7;P- 



de de 



1 d 2 F cx _ 1 d dp cx (e) 
2 pN de 1 de 2 ~ 2 P d(e) de 



(31a) 
(31b) 
(31c) 



with all formulae evaluated in the monodisperse reference system. To clarify the meanings of the various derivatives, 
let us read Eq. (|3l|a) from left to right. It states that the parameter A, which is a property of a monodisperse phase, 
is the rate of change of excess free energy density as the mean property (e) of all particles in the phase is changed. 
The next equality asserts that A is the overall number density p times the rate at which the extensive excess free 
energy changes when the property e of particle number 1 alone is changed. Finally from Eq. (|3l|a), this is the same 
as p times the variation in excess chemical potential for different species. (Though different species are absent from 
the monodisperse system, their excess chemical potentials remain well defined. 
In sections 



IV A 



and IV B two methods are developed, using the various relations in Eqs. (pl|) to ascertain the 



values of A, B and C in any given system. 
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A. Thermodynamic perturbation theory 



In this section, polydispersity is treated as a perturbation to the Hamiltonian of a system. The Hamiltonian of a 
polydisperse system (neglecting the trivial kinetic part) is H(T), a function of the set T of particle coordinates (both 
positional and internal) which defines a configuration of the system. Given that a monodisperse reference system in 
the same configuration T has Hamiltonian H m (T), we define the perturbation H(T) by 

H(T) = H m (T) + H(T). (32) 
By Boltzmann and Gibbs, the Hclmholtz free energy of the polydisperse system is (with units fc^T = 1) 

F = N Jde p{e)(hx[Np{e)] - l) - In J dT exp[-H{T)] 

and that of the monodisperse system is 

F m = N[lnN - 1] - In J dT exp[-H m (r)]. 

In these two equations we may set H and H m to zero to find expressions for the ideal part, which we subtract, leaving 
the excess part of the free energies, 



and 



F cx = - hxj dr exp[-H (T)] + In J dT 
F™ = - In JdT exp[-ff m (r)] + In J dT. 



By applying Eq. (|32|), we see that the difference between these expressions is 

F ox - F° x = - In 



JdTe- H ^e- H 



JdTe- H ™ J 

where the argument of the logarithm is the thermal average of exp(— H(T)) in the monodisperse system. Let us denote 
the thermal average of a stochastic variable C by ((C)). This is a Boltzmann- weighted average over configurations, 
denoted by double brackets to distinguish it from averages over the distribution of species p(e), for which we have 
been using single brackets (. . .). Hence, the excess free energy of a polydisperse system, in terms of a perturbation H 
to a monodisperse system, is 



F cx = F° x - In 



-H 



(33) 



where the subscript m indicates that the thermal average is performed in the monodisperse system. Equation (|33j) 
has the familiar form of thermodynamic perturbation theory but, for polydispersity, it applies only to the excess part 
of the free energy. The ideal parts of the mono- and polydisperse free energies differ by a non-perturbative amount. 

Note that we are not making any approximation of the interactions, unlike many perturbation theories which use 
an ideal or harmonic reference state. Hence, as with our previous analysis, Eq. ( |33] ) treats the effects arising purely 
due to polydispersity, in a fully interacting system. 

We now apply the second equality in Eq. (pip,) to the perturbative expression for the excess free energy, Eq. (|33D , to 
find A by varying the properties of a single particle in an otherwise monodisperse system. The first term in Eq 
-FIS* is a constant and therefore does not contribute. We find 



A/p 



dF c 



_d_ 

dei 



ln((e 



-H 



yielding 



A, P Jim 



\W// 



(34a) 
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where the 'hat' ( " ) has been dropped from the Hamiltonian H since, by definition, the perturbation H is the only 
part that varies with ei. Similarly, for Eqs. (pip) and (pl|c), 



(34b) 



P \9ei9ei// m \\ da // m \\ de x // m \deidei// 

2C _ II 8 2 H \\ //9H\ _ //MLdH\ 

P N \ de,de 2 // m + \\ de x t\ de 2 /) \ de t de 2 // 



(34c) 



In principle, the derivatives of the Hamiltonian appearing in Eqs. (34a), ( p4q ) and ( p4q ) are known from the 
microscopic physics of a given system. For instance, in a system of particles with central, symmetric, pairwise 
additive interactions, whose inter-particle potential is U(r,ei,€j), the Hamiltonian is 



N N 



H = lj2J2u(\r l -r J \,e l ,e J 



(35) 



i=i j=i 



In terms of the matrices of derivatives of the potential, 

8U(r,e,0) 



Ui(r) 



Uu(r) 



U 12 (r) 



de 

d 2 U(r,e,0) 



de de 
d 2 U(r,ex,e 2 ) 



de\de 2 



e=Q 
e=0 

e 1= e 2 =0 



this yields 



A = p z Ui (r) g m (r) Airr 2 dr. 
lo 



(36a) 
(36b) 
(36c) 

(37a) 



Here, g m (r) — ((p(O)p(r))) m / p 2 is the radial distribution function in the monodisperse reference phase. To obtain 
Eq. fl37|a) , we have used the fact that the density of particle centres in the system at any instant is p(r) = ^ (r — 
Ti). Equations (341, 34c) can similarly be evaluated for central, symmetric, pairwise additive interactions, giving 



B = Y J Q [C/ll(r) ~ C/ iW^iW]5m(r)4^r 2 dr + F — J d d rdV U 1 {r)U x {r f ) [g m (r) g m (r') - g£>(r,r') (37b) 
C =yI°° [Ul2{r) - Wi(r)]ffm(r)47rr a dr - ^ J d 3 rdV U^U^g^ (r, r') 



d 3 rd 3 / dV Ui(r)Ui(\r' - r 1 ' 



,(r)g m (\r'-r"\)-g£\r,r>,r") 



(37c) 



where ^ 3) (r,r') = ({p(0)p(r)p(r'))) m / p 3 and gg> (r, r' , r") = ((p{0)p{r)p{r')p(r"))) m / p 4 . Notice that Eqs. (0a, b, c) 
require a knowledge of spatial correlations in the monodisperse system only. In the appendix, it is shown how general 
thermal averages are perturbed by polydispersity, for such systems with soft, pairwise interactions. 

Equations ([m]) and (|37| ) allow the parameters affecting phase equilibria (to lowest order in the width of the distri- 
bution) to be calculated for any system with soft (i.e. diffcrcntiable) interactions. We shall put these expressions to 
use on some real systems in sections [v] and VI . 



B. Scalable systems: The special case of hard spheres 



In this section we shall consider a system of prevalent theoretical interest, polydisperse hard spheres. Hard spheres 
interact via a potential which is zero, except for configurations where the spheres overlap, which are forbidden and 
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hence have infinite potential energy. Hard-sphere systems attract interest because they exemplify substances with 
short-range repulsive interactions, while lacking a characteristic energy scale. This leads to temperature-independent 
behaviour which is purely entropy-driven, and engenders simplicity due to the small number of tunable param eters 



Unfortunately, hard spheres are an exceptional case to which the thermodynamic perturbation theory in section IV A 
cannot be applied, since the internal energy is zero and the interaction potential is discontinuous. Hence the values 
of the parameters A, B and C cannot be found by that method, though their formulations in Eqs. (|3l|) still hold. 

In this section we shall exploit a special property of the hard-sphere system to establish the exact value of the 
parameter A (previously approximated J35|). In fact the method can be applied to any system of particles whose 
interactions are 'scalable' — a term which will be elucidated below. With a knowledge of A, the partitioning of hard 



spheres between slightly-polydisperse phases is fully determined in Eq. (|14|), (|15| ) or (18). Unfortunately it is not clear 
how such a method might be used to determine B, the other parameter on which the phase boundaries depend. In 
section 



VI its value is taken from an approximate expression for the free energy of polydisperse hard spheres, of which 
there are many examples in the literature p6|p^j23] | . 

To calculate A, the coefficient of (e) in the free energy expansion (Eq. (g)), for polydisperse hard spheres, we use 
the first equality of Eq. (|3l]a). Let the radius r of each sphere be measured relative to a reference length ro thus: 
r = (1 +e)fo- As this is the only property that varies from sphere to sphere, the value e which characterises the 



particles is a scalar, so Eq. (31a) reduces to a scalar equation. That equation states that A is the rate at which 
the free energy density changes when the first moment (e) of size deviations is increased, while the other moments 
remain stationary^. When this hypothetical change is made, the particles all grow by the same amount, so the system 
remains monodisperse. In terms of the unique radius r of the monodisperse particles, the derivative is 



dF ex dT 



'CX 



'n- 



d (e) dr 

which will be evaluated at r — ro- One can imagine making this change in two stages: firstly the whole system is 
scaled up, so that the particles, and the space between them, and the volume of the system all increase, while the 
concentration 4> (the fraction of space occupied by particles) is held constant. Then the system is compressed back 
to its original volume while the particles remain at their new large size, so the concentration increases. The resulting 
change in the extensive excess free energy (F cx ~ VJ- CX ) is given by 

dF cx \ fdF cx \ fdF cx \ ( dV 



dr Jv,n \ dr J<t>,N V dV J r N \dr /cj)N 

The last term here is due to the reversible work done against excess pressure when the system is compressed, with 
{dV / dr)^^ — 3V/r being a conversion factor between length and volume. So we have 

' ■ > ' > +3 . 



\ dr J V ^ N \ dr 7 ^ N 

The first term on the right hand side is trivial to calculate because there exist no length scales (such as the range of 
an interaction) that remain fixed as the particles grow. Hence the growth is simply an overall scaling jL7|. This is 
what was meant by the term "scalable system" above. By writing the free energy in terms of the usual configurational 
integral over particle positions measured in units of the thermal de Broglie wavelength, it is easy to establish that 
this term vanishes, since the excess part does not vary with the overall scale factor. Hence, for hard spheres p7y , 

A = 3P CX = 3(P-p) (38) 

where, as usual, ksT has been set to unity. Since, at equilibrium, P is a constant for all phases, the coefficient in 
Eqs. p5| , |l8| ) may be written [A/p]^ = 3P[p -1 ]^ for size-polydispersity in scalable systems. 

Though it is not apparent that an analogous method exists for calculating B, the parameter C is calculable from 
the first equality of Eq. ([Il]c), by taking the second derivative of excess free energy density with respect to particle 
size in a monodisperse scalable system. We shall not perform that calculation here. 



5 In a monodisperse system, the distribution is a delta function. It is therefore not possible to vary the mean while holding 
the other moments constant. Nevertheless, being of higher order in the small quantity e, the higher moments are stationary, 
so the partial derivative is well defined. 
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For monodisperse hard spheres, there is a transition |38],[48| from a crystal of volume fraction 0.545 to a fluid 
phase of volume fraction 0.494 which, according to the Carnahan-Starling equation of state |3{J, exists at a pressure 
6.17 ± 0.02, in units scaled by fc^T and by the volume of a particle. Hence, from Eq. j38|), the coefficient [A/ p]^^ a , 
which appears in Eq. (18), evaluates to 3.51 ± 0.04. This is consistent with simulation results in which the average 
sizes of polydisperse hard spheres were measured in coexisting fluid and crystal phases, and the difference plotted 
against the variance of the overall distribution ||(J. The gradient was found to be 3.55 ± 0.01, and changed by only 
a few percent up to polydispersities of 6% and more. 



V. EXAMPLE: CHEMICALLY POLYDISPERSE POLYMERS 



It has been necessary to plough through a fair amount of mathematics in order to arrive at some simple and practical 
formulae. The convenience of this approach is that the lengthy derivations have been performed once and for all. 
They will not require repetition for each application. Lest the reader has lost sight of the wood for the trees, let us 
demonstrate how easily the formulae may be applied to a model of polymeric phase equilibria. 

We consider a Flory-Huggins-style pl| description of a polymer blend. The reference system is a binary mixture of 
polymers (species l A' and '£?') of equal size L, but two different chemical types. The chemical nature of each species 
is parameterised by a number in the interval (—1, 1), which may be interpreted as the hydrophobicity of the molecule. 
We make species B chemically 'neutral' (hydrophobicity zero), while the value for species A is a. If the concentration 
of .A-polymers is </> then that of £>-polymers is 1 — 4>, and the mean-field free energy density is given by 

LT m = <j> In <f> + (1 - </>) ln(l - <j>) - X <t>\ (39) 

where <j)\ = a<p is the 'hydrophobic concentration'. The Flory-Huggins interaction parameter \ determines the 
strength of attraction between chemically similar species. This model system separates into coexisting „4-rich and 
„4-poor phases for \ = %a 2 > 2, the concentrations of which are easily found e.g. from a double-tangent construction 
Jil"| on the free energy (Eq. (39)). The resulting phase diagram in the (4>, x)-plane is shown by the solid line in Fig. [I]. 



If component A is now made chemically polydisperse, by varying slightly the constituent monomers on each polymer 
molecule, while component B remains monodisperse and chemically neutral, acting only as a polymeric solvent, then 
the overall hydrophobic concentration of polymer A is the mean of a distribution: 



4>i = <fi I ap(a) da — (f> (a) 

and the free energy density becomes 

LT = (1 - 4>) ln(l -</>)- x fi + I <t>p(a) ln(<£p(o)) da - ( 40 ) 



The model free energy of Eq. (40) was studied previously to demonstrate a sophisticated approximation scheme 
for systematically reducing the dimensionality of polydisperse phase equilibria problems [^3| . Accurate cloud (o) and 
shadow (*) points were calculated by that method, and are reproduced in Fig. [j] (where X = X a o) f° r a Gaussian 
parent with (a) p = ao = 0.5 and polydispersity a — 8% fl^ . 

For comparison, our lowest-order perturbative formulae for the cloud and shadow points (Eqs. ( p6| , p7f)) can be 
evaluated with ease. Note that the last term of Eq. ( |40| ) has the form of an ideal free energy density if we identify 
cj) with p in Eq. (Q) (up to an irrelevant term linear hip). Hence, by writing a — ao(l + e) to expand about the 
monodisperse reference system with a = ao = 0.5, Eq. ( f40| ) can be cast in the form of Eq. (||), with 

A = -2x0 2 
B = 

c = -x4> 2 

all of which are scalar quantities since only one property (hydrophobicity) is polydisperse. The cloud and shadow 



points (<j) c — 0™ + 5 C , cj) s — + S s ) are obtained by substituting these values into Eqs. (26, |27J) together with the 
isothermal compressibility derived from Eq. (|3S|), «;=(! — 0)/(</> — 2x4> 2 {l — </>)), yielding 



2a 2 Y i£Js Z£_^Z£ ZjlJZ±. (41) 



2-2 C(i-0(C-0(2^ 



2a z y- z s s c s (42) 
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where ffi 1 , are the respective points on the monodisperse binodal (solid line in Fig. [j]). 

Eqs. (|l|, |§) are plotted in Fig. as dashed and dash-dotted lines respectively. As expected, the perturbative 
expansion scheme breaks down near the critical point, where the compressibility diverges and the separate reference 
phases disappear. This invalid regime extends over only a small range of x, and the phase boundaries rapidly converge, 
with increasing \, to a very accurate answer. 

Notice that the miscibility gap is broadened by polydispersity. Equation ( |29| ) shows that the polydisperse pertur- 
bation to the gap must be positive whenever B = B' = 0, i.e. when there is no explicit dependence of the excess free 
energy on the variance of the distribution of species. 






0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 

FIG. 1. Phase diagram for the chemically-polydisperse Flory-Huggins model with a chemically neutral polymeric solvent, in 
the (0, x)-plane. Solid line: monodisperse binodal. Dashed line: perturbative cloud curve (8% polydispersity). Dash-dotted 
line: perturbative shadow curve. Circles: exact cloud point. Stars: exact shadow point iflif . 



We also have formulae prescribing the chemical properties of the two phases. By definition, the distribution of 
polymers in the cloud phase is just the same as the parent distribution that was input to the system. The shadow 
phase, on the other hand, has chosen its own preferred mix. The mean hydrophobicity of type-,4 polymers in this 
phase is given by Eq. (|l|a) (or equivalently Eq. @) as (a) a = a (l + (e) s ) = a (l + 2cr 2 x(</> s - <t> c ) + 0((T 3 )), which 
is also in good agreement with exact numerical data. 



VI. EXAMPLE: STERIC ALLY-STABILISED COLLOID-POLYMER MIXTURES 



To further illustrate the utility of the formulae derived in this article, let us make a case study of a familiar complex 
fluid: a colloid-polymer mixture p4j| . The phase behaviour of the monodisperse system is first briefly reviewed. A 
theoretical understanding of its behaviour is well-established Jig] . However, the effects of colloidal polydispersity on 
the system have not previously been calculated, despite the fact that the results of colloidal synthesis are inevitably 
slightly polydisperse. The effects of polymeric polydispersity on the system, which have been modelled previously 
J46|, will not be analysed. 



A. Established monodisperse behaviour 

In this system, small balls (typically of diameter 10 nm to 1 (im) of PMMA ('Perspex') [Q or silica |47| are 
suspended in an organic liquid. To avoid the balls (or 'colloidal latices') clumping due to van der Waals attraction, 
they are each given a short brush-like coat of 'steric' polymer chains, making the latices behave as hard spheres, with 
no long-range interactions but a very short-range strong repulsion. Hard spheres of this kind are known |48| ] to have 
a simple phase diagram, and carry no latent heat. The only controlling parameter is the fraction of space d> occupied 
by the spheres. A single phase transition exists, between an amorphous fluid state and an FCC crystal, which coexist 
at <f» « 0.494 and d> ~ 0.545 respectively [Q. The phase behaviour is enriched by the addition of free polymer chains 
to the mixture. Each of these chains, being very long, bunches itself into the form of a self-avoiding random walk. 
Such a 'ball of string' cannot penetrate a colloidal latex, and so interacts with the colloid particles as if it too were 
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a hard sphere. On the other hand, when two randomly coiled polymer molecules meet, they can interpenetrate and, 
with a suitable choice of suspending solvent, can be made almost entirely non- interacting, like an ideal gas Q,fl5]. 
The ideal gas pressure that they exert on the surface of an isolated colloidal hard sphere is isotropic. However, if a 
second hard sphere is so close to the first that the non-penetrating polymer coils cannot fit between them, then the 
polymeric pressure is not felt on the adjacent parts of the colloids, so they are pushed together. 

Asakura and Oosawa J5(J showed that this mechanism could be described in terms of an effective attraction between 
colloidal spheres, known as the depletion interaction. Thus, rather than considering a mixture of two components 
(colloid and polymer) in the solvent, one can imagine a single-component system, consisting of spheres with a hard 
core repulsion and a longer-range pairwise[] attraction. This effective attractive potential U (r) between two latices of 
radius a and centre-to-centre separation r is simply the product of the (ideal) polymeric osmotic pressure Tl p and the 
volume of space from which the centres of polymer coils are excluded by the hard spheres [Q . At r = 2a (contact of 
the two hard spheres) the potential rises discontinuously to infinity, and for r > 2 (a + r g ), where r g is the radius of 
gyration of a polymer chain, there is no interaction, since the diameter of a polymer coil (2r g ) sets the range of the 
effective potential. 

The phase diagram of this system in the (</), n p )-plane has been calculated for the monodisperse case |j45| , using 
results from the Percus-Yevick integral equation theory for hard spheres fl52| . It was found (in agreement with 
experiments ifhj ) to resemble the phase diagram of simple atomic substances such as argon, with n~ 1 playing the 
role of temperature. As well as the colloidal crystal, there are low- and high-concentration amorphous phases, dubbed 
'colloidal gas and liquid', which meet at a critical point. We shall study the 'fluid-fluid' coexistence of these amorphous 
phases, thus avoiding problems of non-ergodicity arising in the polydisperse crystal [pi]. 



B. Calculating the polydisperse phase equilibria 

The strength of the treatment of polydispersity in this article is that it uses, as input, properties of the monodisperse 
system. The phase behaviour of monodisperse colloid-polymer mixtures is well understood, so application of the new 
formulae for polydisperse phase equilibria will be straightforward. In addition to the monodisperse phase diagram, 
we require the functions A(p) and B(p), defined by Eq. (||), which are scalars in this case since the particle radii alone 
are polydisperse. 

Unlike the previous example in section we are not given an expression for the free energy of the polydisperse 



system. According to Eq. (J31 
given in Ref. E{3]. It is simp' 



a), A(p) can be found from the monodisperse free energy, for which we use the expression 
y the derivative of the free energy with respect to the radius of monodisperse colloidal 
particles (scaled by the radius), at fixed number density. Hence an expression for A(p) is straightforwardly found, so 
long as care is taken not to overlook the dependence of concentration on particle size, <f> = ^ira 3 p, and to keep the 
polymer size r g fixed while taking the derivative. 

It is not so trivial to evaluate B(p) without a prescribed polydisperse free energy. Firstly, we must specify the 
difference between the Hamiltonian of the system with slightly polydisperse colloid, and that of a monodisperse 
reference system. Let us write the Hamiltonian of the colloid-polymer mixture with polydisperse colloid as a sum of 
two parts: 

H = Hft S + fff p (43) 

the Hamiltonian i/jjs of the hard-sphere colloid, plus effective interactions due to the depletion of polymer. 

Following the analysis p5] of the monodisperse case, we assume that the resulting free energy is also separable into 
two parts: 

F = F HS + F dcp (44) 

where Fgs is approximately the free energy of polydisperse hard spheres in the absence of polymer. This is consistent 
with the Weeks-Chandler- Andersen (WCA) approximation Q which states that the spatial distribution of repulsive 
particles (in this case hard spheres) is not significantly altered by the introduction of additional interactions which 



In particular, the second virial coefficient of the polymeric osmotic pressure can be made to vanish, though higher virial 
coefficients remain p9| . 

7 A more accurate description would include corrections to the pairwise attraction, since three or four latices in close proximity 
exclude polymer coils from a region with non-trivial geometry. 
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are purely attractive (here contained in ). Accordingly, the polymeric contribution to the free energy -Fd ep will 
be approximated by averaging the depletion Hamiltonian -ff|ep over the pure hard sphere distribution. 
Note that the structure of Eq. (Q) implies that the desired function B(p) also splits into two terms, 

B = Bus + B dcp (45) 



The depletion interaction is soft, so the thermodynamic perturbation theory of section [V A is suitable for calculating 
its contribution in -Bdcp- The hard sphere potential, on the other hand, is non-differentiable and therefore inappropriate 
for description in terms of an energetic perturbation, as noted above. Other methods will be used to deal with the 
hard sphere part Bus- 

Let us tackle the polymeric (depletion) contribution first. As stated above, the effective Hamiltonian (i.e. the 
polymeric free energy for a fixed configuration of colloidal spheres) H^ p is simply the product of the ideal polymeric 



osmotic pressure and the volume of space from which centres of polymer coils are excluded. Let us define the polymer- 
s/a, and v s = | 



colloid size ratio as £ = r g /a, and v s = §7ra 3 to be the volume of a reference sphere with the mean radius. An isolated 



colloidal sphere of radius a(l + e) excludes polymer coils from a volume v s (l + £ + e) 3 . Hence, even for well-separated 
colloids beyond the interaction range, there is a bulk contribution to the depletion Hamiltonian, 



N 



#buik = ^^(l+e + e,)" (46) 
^ i=i 

where the polymeric osmotic pressure has been re-scaled: Ii p = v s ^ 3 H p , so that the (</>, n p )-plane will correspond to 
phase diagrams presented in Rcf. for colloid-polymer mixtures. For monodisperse colloid, -ffbuik is usually ignored, 
since it is a constant, independent of volume fraction (j>. For a polydisperse system, on the other hand, it cannot be 



neglected, as it contains e-dependence. From Eq. (34b), its contribution to B(p) is 



Bbuik = 3fV(l + £)/e 3 - (47) 

When two colloidal latices, of radii (1 + e\)a and (1 + £2)0, have a centre-centre separation (2 + t\ + £2)0 < r < 
(2 + 2£ + t\ + £2)0, their individual regions of polymeric exclusion overlap, so that the total volume available to 
polymer increases, giving rise to the depletion attraction. It can be shown by some elementary geometry that the 
overlap volume of the two spheres of depletion; radii n = (1 + £ + ei)a and T2 — (1 + £ + £2)0; is 



Vr 



overlap 7T 



(rl+r%)r 2(r\ + rf) (r? - ^ 2 



'27 



12 2 3 4r 



The resulting interaction potential U(r, e%, £2) = — n p V over i ap gives a pairwise contribution to the effective Hamiltonian, 
as in Eq. (Bq). Its derivatives, as defined in Eqs. (|3q), are 



U i M = " ^TTT^ ( 2 + 2 ^ - r M n p ( 48a ) 



3(1 + 0. 

C/iiW = ^(^-4[l + e] + 2^[l + C] 2 ) U P (48b) 
which can be substituted into Eq. (p?p) for -Bdcp- Additionally, that equation requires the pair and three-point 

(31 

distribution functions g m (r), gin (r,r ) for the monodisperse system. At the level of the WCA approximation, these 
are replaced by the pure hard-sphere (HS) distribution functions, g « ff HS , for which we use the Percus-Yevick 



expression ]54j]. As a simplifying assumption, let us also use the Kirkwood superposition approximation 55 for the 
three-body correlations, 



,(3) 



(r,r') w g(r)g{r')g(\r - r'\). 



It remains only to find the hard-sphere contribution to B. Unfortunately, the methods developed here do not enable 
a first-principles derivation of this quantity, so we must look to other analyses of polydisperse hard spheres for an 
evaluation of -Bjjs- In particular, we refer to the free energy expression for the polydisperse hard-sphere fluid due to 
Boublik, Mansoori, Carnahan, Starling and Lealand J36| (BMCSL), which is known to reduce to the Carnahan-Starling 
free energy J3l| in the monodisperse limit. The application of this expression perhaps requires some clarification. The 
use of results from other studies of polydispersity by no means depreciates the present analysis. In cases where a 
polydisperse free energy is already known, one is spared the application of first-principles methods such as section 
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IV. However, the results of section III (or equivalent methods) are still required, in order to derive (to lowest-order) 
the phase equilibria from that free energy. Indeed, the gap between a knowledge of the free energy and a knowledge 
of the phase equilibria is evidenced by the BMCSL free energy, whose regime of thermodynamic instability was only 
recently established |20| . 

To extract the required function, the BMCSL free energy is cast in the form of Eq. (|^) by expanding to order e 2 , 
at a density p — 4>a/v s , giving 



T7CX _ 



^g(4-3fo) 
(l-0o) 2 



(£) 6^(2-00) +(e2)3 ^ o 



(1-0O) 3 

fo(l + 2fo)(3 + <ft)-^) 
(1 - </>o) 4 



00(1 + 30o -2$) 



ln(l - fa) 



+ ln(l - 0q) 



0(e 3 ). 



(49) 



The zeroth-order term is the excess free energy of a monodisperse hard-sphere fluid, and recovers the Carnahan- 
Starling equation of state. By comparison with Eq. (||), the coefficient of (e 2 ) is the required function v s Bns- 

For definiteness, let us consider a polymer-colloid size ratio £ = 0.4. For a given polymeric osmotic pressure 
tl p , the coexisting concentrations <fi of monodisperse colloidal fluid phases are given in Ref. EEJ. These values for 
the concentrations of coexisting monodisperse gas (</> m ) and liquid (4> l m ) phases are substituted into the expressions 
calculated above to obtain values of A and B for the coexisting phases. The monodisperse gas-liquid phase boundaries 
0m) 4>m are shown as a solid line in Fig. ^ from the critical point at n p w 0.41 to the triple point at Tl p rj 0.54. 

The calculated values of A and B and their derivatives are used in Eqs. (|2(], [2?]) to find the change in each phase's 
cloud- and shadow-point densities due to polydispersity, which are translated into concentrations via Eq. (^8|) . These 
shifts in concentration are added to the monodisperse phase boundary in Fig. || to give the resulting cloud and shadow 
points C,S = m s + S(j) c ' s for polydispersity a = 8%, shown as dashed and dash-dotted curves respectively. We see 
from the figure that the expansion has remained well controlled (has not blown up) even for values of n p very close 
to criticality of the reference system. The theory predicts a widening of the coexistence region (the miscibility gap) 
so that liquid will condense from a polydisperse gas of lower concentration than in the monodisperse case. We also 
observe that when the gas phase exists in an infinitesimal quantity (at the shadow point), it is more concentrated 
than at its cloud point. The same applies to the liquid phases. Similar results are also found for different size ratios 
£ (studied for 0.3 < £ < 1). 




I i I i I i I i I i I 

0.1 0.2 0.3 0.4 

colloidal concentration 

FIG. 2. Colloidal gas-liquid phase boundaries for colloid-polymer mixture with size ratio £ = 0.4, in the plane of colloidal vol- 
ume fraction <f> and polymeric osmotic pressure lip. Solid line: monodisperse binodal. Dashed/dash-dotted lines: perturbative 
cloud/shadow curves for 8% polydisperse colloid. 



The above calculation of B(p) for the colloid-polymer mixture involves a degree of inexactitude, since the radial 
distribution function is approximated by the Percus-Yevick hard-sphere distribution function, while three-point cor- 
relations are even more crudely estimated. To assess the importance of precision in the correlations, the cloud and 
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shadow points were re-calculated using a Heaviside step function g(r) = 0(r — 2a) in place of the Percus-Yevick radial 
distribution. The polydispersity-induced shift in the high-density cloud point changes sign when the correlations are 
neglected, indicating, not surprisingly, that the correct correlations are important in determining the high-density 
phase boundary. Though the positions of the other phase boundaries are also affected by the change in g(r), certain 
qualitative features remain unaltered, and are therefore expected to be reliable predictions of our approximate model. 
In particular, each shadow point remains more concentrated than the equivalent cloud point and, far from the critical 
point, the gases are shifted to lower concentrations by polydispersity. 

According to Eqs. and (|l8|), the difference between the normalised populations of the coexisting phases is 
proportional to — [A/pJ^and the fractional difference in mean particle size is [(e)]g = — pf g o\ where a p is the 
overall polydispersity. So the coefficient — [A/p]^ determines the amount of partitioning or fractionation. It is plotted 

for the gas-liquid phase boundary (from tl p = 0.41 to 0.54) in Fig. |j. Note that it is positive, indicating that, on 
average, the liquid phase contains larger particles than the gas. The hard-sphere contribution is negative, but is 
outweighed by the depletion part. 
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polymeric osmotic pressure 

FIG. 3. Fractionation per unit variance ([(e)]^ / &p = — [A/p]g) against polymeric osmotic pressure Ii p for gas-liquid coexis- 
tence of a colloid-polymer mixture with size ratio £ = 0.4 and colloidal polydispersity a p — 6%. 



VII. SUMMARY 

The results derived here are exact in the limit of low polydispersity, and can therefore provide accurate informa- 
tion on a great many systems, whose manufacture attempts to approximate the ideal of monodisperse constituents. 
Additionally, the results will be of pragmatic use even in systems with a wide scatter of particle properties, in pro- 
viding qualitative predictions and estimates for the effects of polydispersity. In any case, minimal effort is required 
to substitute values into the formulae. 

Equations ( [l4| ) and, equivalently, (|l5|) show that, in the narrow limit, the difference (due to fractionation) between 
normalised distributions in coexisting phases does not depend on the volumes of the phases, and is in fact remarkably 
simple, requiring only one system-dependent parameter per polydisperse property. From Eq. (|l8|), the mth moment 
about the centre of the distribution differs, between coexisting phases, by an amount proportional to the (m + l)th 
moment of the parent^. For hard spheres, the constant of proportionality has been determined (Eq. (|38|)), and agrees 
with data from simulation p0| . The calculation is indistinguishable from the data at 2% polydispersity, and deviates 
by only 5% at 4% polydispersity. 



plus terms of the order of the (m + 2)th moment, which become important in a near-symmetric distribution if m is even (BU 
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Additionally, at the nematic cloud point of the polydisperse Zwanzig model of hard rods |?3| , the form of Ecj 



has been shown [|56| to hold to within 5% up to the remarkably large polydispersity of 50%. In section |III A| 7 
was shown that this simple rule can lead to 'convective fractionation' for multiply-polydisperse systems, whereby one 
property is partitioned between phases due to a driving force on another. 

The expressions found for the shift in the phase boundaries due to polydispersity, at a general coexistence (Eq. (|2 



and at the cloud (Eq. fl2q)) and shadow points (Eq. (|27[)), show that the shift in the density of a phase is proportional 
to its isothermal compressibility and to the square of the polydispersity. This explains why these phase boundaries 
are parabolic near a = when plotted as a function of polydispersity a in a range of models; see e.g. |]6|,[7|, p3|j24|j3^ , |3^| . 
The other parameters on which these shifts depend can be extracted directly from the form of the polydisperse free 



energy if it is known or, if not, from the Hamiltonian via the thermodynamic perturbation theory in section IV A. 
In the case of pairwise interactions, the relevant parameters are obtained from a knowledge of only two- and three- 
body correlations in the reference system, even when larger clusters are responsible for the existence of the phase 
transition in question. This illustrates a strength of the small-polydispersity expansion: that the underlying reference 
system takes care of complicated many-body interactions, so they do not need to be calculated in the analysis of the 
polydispersity. 



The application of the methods to the fluid-fluid coexistence of colloid-polymer mixtures in section VI revealed that 
colloidal polydispersity leads to a widening of the coexistence region, and favours larger particles living in the liquid 
phase. Similar methods are used in appendix |A| to find how polydispersity alters correlation functions. 

Other methods exist for calculating polydisperse phase equilibria. Whatever the theoretical formalism by which 
one chooses to analyse a polydisperse system, the results, if correct, will tend in the limit of low polydispersity to 
Eq. ( p"5| ) for the degree of fractionation, and Eq. ( psj ) for the movement of the binodals. 
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APPENDIX A: THERMAL AVERAGES AND SPATIAL CORRELATIONS 

In polydisperse systems, correlation functions are easier to calculate than phase boundaries, as there is no parti- 
tioning of the sample to consider. The thermal average of a stochastic quantity ( in a single polydisperse phase is 
given by the usual perturbative expression with respect to a monodisperse reference system 

« C )) = «CO)m (A1) 

«e-"» m 

We shall require the second-order expansion of this expression, in the small quantity H , which must therefore be small 
(in units of ksT), 

((C)) = «C»m - (l + {(H)) m ) («C#» m - ((C» m «#»m) + \({{CH 2 )) m ((0) m ((H 2 » m ) + 0(H 2 ). (A2) 
Let us restrict the discussion to a system of particles interacting via a pairwise-additive, symmetric, isotropic, central 



potential U(r, £1,62), for which the Hamiltonian is given in Eq. (35). With the Hamiltonian Taylor-expanded to 
second order in e, the perturbation can be written 



N N 

i=l j=l 

in terms of the derivatives of the interaction potential defined in Eqs. (|36|). Substituting this expression into Eq. (|A2|), 
we keep terms to second order in e, and measure the deviations e with respect to the mean properties of particles in 
the single-phase system, so that (e) = 0. This yields 
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((C)) = «C»„ - (ee) :^E{((C *(k< ~ r iD»m - ((0) ro - r,|))) m } 

N N N 

+ («) : 2EEE{«^d r '- - ^D» m - «0) m - r-iDC/idn - r fc |))) m } (A3) 



i=l j=l k=l 



where 



*(r) = I7u(r) - iz7!(r)t/i(r) 



The thermal averages of functions of the stochastic particle positions can be re- written in Eq. (A3) in terms of the 
number density field in the monodisperse system, which at any instant is p(r) = ^ 3 H r — r i); where 5^ is the 

three-dimensional Dirac delta function. Finally, the thermal average of a given stochastic variable £ becomes 

((C)) = «C» m - ~ (ee) : J dr' dr" ( ({( p(r')p(r"))) m ((C)) ra ((p(r')p(r"))) m ) *(|r' - r"\) 

+ \(ee): J dr' dr" dr'" ( ((Cp(r>(r")p(r'"))) m - «C)) m (WW'W"))) m ) U x (\r' - r"\) Oi(|r' - r"'\). (A4) 

A case of interest, for instance, is C = p(0)p(r). Then Eq. (A4) gives the perturbation, due to polydispersity, of 
the pair correlation function for total density as 

«p(0)p(r)» - «p(0)p(r)» m - ~ {ee) : / dr'dr" f ({p(0)p(r)p(r')p(r"))) m - ((p(0)p(r)» m ({p(r')p(r"))) m ) *(|r' - r"\) 



+ \ (ee) :|dr'dr"dr'" ( ((p(0)p(r)p(r')p(r")p(r"'))) m - ((p(0)p(r))) m {{p(r')p(r")p(r"'))) m ) C/xQr' - r"|) Ox(|r' - r"'\) 

(A5) 

which depends, to second order in the standard deviation of species, on four- and five-point correlations in the 
monodisperse reference phase. 
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